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Abstract. - We present a model for the dynamics in energy space of multicanonical simulation methods 
that lends itself to a rather complete analytic characterization. The dynamics is completely determined by 
the density of states. In the ± J 2D spin glass the transitions between the ground state level and the first 
excited one control the long time dynamics. We are able to calculate the distribution of tunneling times 
and relate it to the equilibration time of a starting probability distribution. In this model, and possibly in 
any model in which entering and exiting regions with low density of states are the slowest processes in the 
simulations, tunneling time can be much larger (by a factor of O(N)) than the equilibration time of the 
probability distribution. We find that these features also hold for the energy projection of single spin flip 
dynamics. 



The simulation of models with complicated energy landscapes, such as spin glasses or models of 
protein folding, has always proved impractical for conventional canonical Monte Carlo methods [1]. 

In a Monte Carlo simulation of a system of N degrees of freedom, at fixed temperature T, energy 
fluctuations about a mean energy E(T) are of order fcsTV 'Nc v . As a result, only states in a range 
of energies (E) — 5E < E < (E) + SE with 5E ~ IcbT\/Nc v are accessible in the simulation. 
In systems with complex energy landscapes, phase space in this narrow energy range breaks up into 
many regions, connected only by states requiring energy fluctuations AE ^> 8E; tunneling times 
between these regions become too long to retain any hope of achieving the asymptotic distribution in 
a reasonable simulation time. 

To overcome this problem, one would like to broaden the energy range of the states sampled in 
a simulation, forsaking the canonical ensemble at a fixed temperature. A variety of methods have 
been proposed to implement this idea, such as entropic sampling [2], multicanonical Monte Carlo 
[3,4], simulated and parallel tempering [5-7], Wang-Landau sampling [8-10], broad histogram and 
transition matrix methods [1 1-14]. In some cases, like multicanonical Monte Carlo and Wang-Landau 
sampling, the aim is to sample all energy levels of the spectrum with equal probability, producing 
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flat histograms in energy space. While this allows overcoming large energy barriers in systems with 
rugged energy landscapes, critical slowing down can remain a problem [14, 15]. 

Given the practical impossibility of measuring equilibration time of an initial distribution, the 
performance of these methods is generally assessed by measuring the so called tunneling time [15-17], 
the time required to cross the entire energy spectrum, from ground states, to states with the energy of 
maximum density of states, and back, or from ground state to anti-ground states, states of maximum 
energy. 

In this letter we focus on the dynamics of this random walk in energy space. We propose a modi- 
fication of the directed network of non-zero transition probabilities which underlies the Markov chain 
of a Monte Carlo simulation, which, while retaining the long time behaviour of the energy projection 
of conventional simulations, lends itself to a very complete analytic description. Our three most im- 
portant results concern the ± J 2D spin glass: (a) the transitions between ground state level and first 
excited one completely control the long time dynamics of simulations in this model; (b) the equilibra- 
tion time of a simulation can be considerably shorter than tunneling time; (c) these conclusions also 
hold for the case of conventional single spin flip (SSF) dynamics. 

To calculate thermodynamic averages in any system, we need only specify the set of points of 
phase space and the corresponding probabilities. To solve this problem using Monte Carlo, we impose 
on this structure a directed network of non-zero transition rates which defines the Markov chain used 
to generate the asymptotic distribution. The traditional SSF connects each point of phase space to 
N first neighbours. The SSF algorithm is, usually, local in energy, i.e., energy differences between 
connected phase space points being of O(l), independent of system size. It has been claimed that 
differences in the scaling behaviour of the ± J 2D spin glass model and a fully frustrated model, both 
with extensive ground-state entropy, reside in the restrictions imposed by the network defined by SSF 
dynamics. Still, the scaling laws of relevant time scales with system size found with our dynamics [16] 
are also similar to those obtained with SSF dynamics [15]. 

The simplest choice that circumvents the problem of local minima is to construct a Markov chain 
that connects each point of phase space to all the remaining points that have the same or adjacent 
energies. In this way, we preserve the locality in energy space of the standard SSF network and avoid 
the need for long paths to connect states that are adjacent in energy. 

The procedure in best explained in reference to a specific model like the ± J 2D Ising spin glass. 
The corresponding energy levels, E i7 and degeneracies, N(Ei), can be calculated exactly using the 
program of Saul and Kardar [18]. Non-zero transition rates connect states with energy Ei to all states 
with energy Ei or Ei±\. The simplest choice of transition rates is to propose the final state uniformly 
from the set of allowed states, and accept the proposed state with a probability which ensures that the 
asymptotic probability of each state equals 1/N(E). Whereas in a SSF simulation most transition 
probabilities to states of the same or nearby energy are zero, in our model they are all non zero, with 
values adjusted to ensure the same equilibrium distribution. This procedure is akin to an averaging of 
transition rates. While this could lead to significant changes in the dynamics, we nevertheless find that 
the most relevant features of the long time SSF dynamics in energy space still hold in our model. The 
advantage of the current model is that the projection of the multicanonical Markov chain in energy 
space remains a Markov process, allowing us the use of a set of analytic tools to study the dynamics. 

We denote by p a (t) the probability of being in a state a; Ni is the degeneracy of energy level Ei 
and Mi = iVj_i + Ni + N i+i — 1 is the number of states accessible in a direct transition from any 
state of energy Ei (the choice of tentative final states excludes the current state). In equilibrium, the 
flat energy histogram condition requires, for a state a of energy Ei, 

Pa = ME7) - N (1) 
For reasons of efficiency, we exclude the negative temperature region: for energies greater than E m , 
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the energy at the maximum of density of states, we replace N(E) by N(E m ) in eq.^ The transition 
rate from a state a to an accessible state (3, iop a , satisfies detailed balance, with the equilibrium distri- 
bution given in eq.Q as usual, it is written as a product of a proposal probability, which we take to be 
uniform among all possible final states, and an acceptance ratio, determined by the detailed balance 
condition: 

1 • (, M(E a )N(E a ) \ 

In the master equation, the transition probabilities u>0 a do not change if we vary j3 (or a) within the 
same energy level. This makes it possible to write the following master equation for the probability of 
having energy E i7 Pi = Y, a PaS Ea ,Ei 

Pi{t+1) = rVPi+iCtJ+IWi-itt) 

+ {l-Ti-Ti^Piit) (3) 
= ^SUjPsit) (4) 



with 



r ' = T7 mm !) TTlTr • ( 5 ) 



U.i ^ V ' M,N, 



This equation defines the projection of the original Markov chain onto to the energy variable; the 
random walk in energy space remains a Markov chain with transition probabilities defined by the 
density of states of the original spin model. This is,of course, an important simplification afforded 
by our choice of transition rates. With the exception of special models, like the infinite range Ising 
model [19], the projection of conventional dynamics in energy space is non-markovian and memory 
effects can have a significant impact on the dynamics [14, 15]. Nevertheless, in the model under study, 
we find below that important aspects of the long-term SSF dynamics are preserved. 
Assume the system starts from an energy level r at time t = 0: 

Pi(0) = 5i, r 

and let Qi(t) now refer to the probability of having energy E{ at time t, given that the system has never 
visited energy level s (s > r). Qi(t) satisfies the master equation in eq.[3]for { < s, with the same 
initial condition as Pi(t), namely Qi(0) = 8i. r , and an additional condition Q s (t) = replacing eq.[3] 
for i = s. The probability of first passage in s becomes 

H sr {r) =r s _ 1 Q s _ 1 (r). (6) 

The master equation can be solved using a normal mode expansion, Qi(t) = ^ a 7 /7A* . The di- 
mension of the transition matrix in energy space, f2, scales linearly with system size, not exponentially 
as in the case of the transition matrix of the Markov process in phase space. The diagonalization of 
SI, becomes a manageable problem allowing the calculation of the eigenvalues, A 7 , the left and right 
eigenvectors, g/and and the coefficients a 7 = gj 9ifi- ^ ^ s a l so possible, of course, to di- 
agonalize the transition matrix for the actual probability distribution Pi (t) using the same formalism, 
and access the decay time of the various eigenmodes of the master equation. The largest finite one (the 
equilibrium distribution, P eq , does not decay) is the equilibration time of Pi(t), r eq . 

The distribution of tunneling time measured in our simulations can be written in the form 



H{t) = Y, H o m (r ~ T')H mQ {r') (7) 
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Figure 1 - Comparison between exact result given by master equation approach and the Monte Carlo results for 
two spin glass samples with L = 6 and L = 12. 



where E m is the energy level corresponding to the maximum of N(E). In Fig. ^ we superpose a 
distribution measured in a Monte Carlo simulation with the rates given by eq. [2]with a calculated one 
for two spin glass samples of linear dimension L = 6, 12; the simulation reproduces the calculated 
distribution very accurately. 

The distribution of tunneling time between any two energy levels decays exponentially with a 
time constant given by T max = — 1/ log(A mQX ), where X max < 1 is the largest eigenvalue of the 
corresponding l~i (eq.|6j. For a given spin glass sample, we denote by t u and Td the longest decay time 
for tunneling Eq — > E rn and E rn — » Eq, respectively; H{t) (eq.[7} will decay exponentially with the 
largest of t u and Td for the given sample. The power law decay seen in [15] appears only when we 
aggregate tunneling times from different realizations of the random interactions ± J. If the distribution 
of the largest decay time, in the ensemble of spin glass samples, has a power law tail, p (r max ) ~ T~" ax 
for r max — > oo, we obtain, H{t) ~ J dxp (x) e~ T > x ~ t x ~ v as t — > oo. 

This asymptotic behaviour is directly related to features of the density of states of the ±J spin 
glass. The inverse transition rate from the ground state level to the first excited level is given by (eq. 
0, 

r -l 1 , ^1+^2 

It was noted in [15] that Ni/Nq has a Frechet probability distribution function, in an ensemble of 
spin glass samples; a similar result occurs for tq (Fig. |2}. This fact completely controls the long time 
dynamics of the Monte Carlo simulations. 

In Fig.|5]we plot r eq , the equilibration time, and t u and Td, the decay time constants for tunneling 
from Eq to E m and from E m to Eq, against To. In the samples with larger To, we observe the following 
relations: 



' c q 



TQ (8) 
Tu W T . (9) 



Note that Td can be two orders of magnitude larger than t u ; tunneling from E m to Eq is much slower 
than from Eq to E m . In fact, as can be seen in the inset of Fig. |3](b), the following relation holds 
asymptotically: 

T d ~N b T u (10) 
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Figure 2 - Comparison between the histogram of tq from 10000 uniformly generated samples of ± J 2D spin 
glass with L — 10 and the maximum likelihood fit to a Frechet distribution function obtained from the original 
data. 



where N b ~ O(N) is the number of energy levels between E and E m (recall that in our simulation 
the energy histogram is limited to this energy range). These results are quite simple consequences of 
the existence of a bottleneck in the dynamics due to transitions between Eq and E\ (large To). 
The equation for Po is (for long time we replace Pi(t + 1) — Pi(t) by dPi/dt) 



dPpjt) 
dt 



-r (P (i)-Pi(t)) 



(ii) 



We denote the tunneling time scale between E\ and E m by T\ ; for samples with very large To, we 
assume to T\. For t 3> t%, we will have for i ^ 



N b -1 N b -l 



P eq (1-PoW) 



where P eq — l/N b is the equilibrium distribution. Replacing in eq.[^we obtain a time constant for 
the decay of Po(t) — P eq given by 



To 



Nh 



N b -1 



To calculate the distribution of tunneling time from Eq to E m , H(r), we use the initial condition 
Pi(t) — Si t Q. Defining Pi(t) — Qi(t) + Ri(t), where Qi(t) is the probability of having energy Ei 



at time t given that the system has not reached E m , Q(t) — 1 Qi(P) i s ^ probability that the 

tunneling time is greater than t, Q(t) — drH(r). For t 3> r-y we have 



Po(t) « Qo(t) 

Pi(t) PB Ri(t) 



for i ^ 



This expresses the fact that the system either has energy Eq, and has not tunneled, or has left the 
ground state and has almost certainly tunneled to E m . Therefore 



dQ(t) 
dt 



dQ Q ^ dP 
~dt dT 
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Figure 3 - (a) Correlation between equilibration time r eg and to for lattice size L = 6 to 20. The dashed line 
shows the asymptotic behaviour where r eq = To. (b) Correlation between both t<j and r u with To. The dashed 
line shows the asymptotic behaviour where r u = To. (c) Correlation between Td/Nt and r«. The dashed line 
shows the asymptotic behaviour where Td /JV b = r u . 

Figure 4 - Correlation between Td/Nt and r u obtained from Monte Carlo simulations of the ± J 2D spin glass 
with SSF dynamics. The dashed line shows the asymptotic behaviour where Td/Nt = t u . 



The decay time of Q(t) (and H(t)) is r„ w ro. 

The asymmetry of r^ and r u should by now be obvious. For tunneling from E m to Eq the initial 
condition is -fj(O) = <5i. m and Q(f) = 532= i Qi(t)- F° r times much larger than ri we expect, now, 

Q t ( t )^-QML. fori ^0 

Since dQ(t)/dt = — ToQi(t) {Q(t) only changes due to transitions between E\ and Sq)> we obtain 

dQ(t) = Tp ( 
d< iV & — 1 

/. e. 



T d = {N b -1)T U W 7VfcT u . 

In simple language, these results can be understood as follows. Tunneling from Eq to E m is 
controlled by the process of exiting the ground state level: the system cannot have tunneled if it is still 
in a ground state. On the other hand, a system can only enter the ground state level if it is in energy 
level Ei ; therefore the tunneling rate from E m to Eq has an extra factor P eq — 1 /Ni, corresponding 
to the probability of having energy E\. 

These results provide a very clear explanation of the correlation between tunneling time and 
Nx/Nq found in [15]. The controlling time scale t = 1 + (Ni + N 2 ) /N is clearly related to 
the ratio Ni/N . This difference in time scales t u and Td should be present in other models; namely, 
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those in which the slowest processes in the simulations involve getting across steep entropy changes. 
The 2D ± J spin glass is an extreme case of this behavior, the dynamics being controlled by the tran- 
sition between the two lowest energy levels. One could ask if this result is an artifact of our simplified 
dynamics. In fact, we verified the same behavior with SSF dynamics. Fig.|4]shows the same type of 
plot as in Fig. [3] with averaged tunneling time (averages taken within a sample) in SSF dynamics: for 
long times the relation of eq. ^|is still verified. Notice that the equilibration time scale of a proba- 
bility distribution, r eq , may be much smaller than tunneling time scale, which, as usually measured, is 
dominated by T4. This can lead to a very pessimistic estimate of the time required to reach equilibrium. 

In summary, modelling the dynamics of a multicanonical simulation in a sort of mean-field way, 
by averaging transition rates to states of the same or nearby energies, we were able to define a Markov 
process in the energy variable, reducing the dimension of the Markov matrix to a manageable size. 
Nevertheless, this procedure preserves the main features of the long time dynamics of a conventional 
simulation. A very complete and physically transparent description of the more salient features of 
the dynamics of multicanonical simulations becomes possible. In particular, we clarified the relation 
between equilibration and tunneling time, and the difference between tunneling away from regions of 
low density of states from tunneling into such regions. 
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